MODELLING AND SIMULATIONS OF MULTI-COMPONENT LIPID 
MEMBRANES AND OPEN MEMBRANES VIA DIFFUSIVE 



Abstract. In this paper, phase field models are developed for multi-component vesicle mem- 
branes with different lipid compositions and membranes with free boundary. These models are used 
to simulate the deformation of membranes under the elastic bending energy and the line tension 
energy with prescribed volume and surface area constraints. By comparing our numerical simula- 
tions with recent experiments, it is demonstrated that the phase field models can capture the rich 
phenomena associated with the membrane transformation, thus it offers great functionality in the 
simulation and modeling of multicomponent membranes. 

1. Introduction. Lipid vesicle membranes are ubiquitous in biological systems. 
Studies of vesicle self assembly and shape transition, including bud formation [24U25I 

and vesicle fission are very important in the understanding of cell functions. 
In recent experimental studies, multi-component vesicles with different lipid molecule 
compositions (and thus phases) have been shown to display even more complex mor- 
phology involving rafts and micro-domains [2] • There are strong evidences suggesting 
that phase segregation and interaction contribute critically to the membrane signaling, 
trafficking and sorting processes 3 . In the literature, the geometric and topological 
structures of multi-component vesicles have been theoretically modeled by minimizing 
an energy with contributions of the bending resistance, that is the elastic bending en- 
ergy, and the line tension at the interface between different components or the phase 
boundary [41 1201 l2lfl I26| . The elastic bending energy first studied by Canham, Evans 
and Helfrich |1 II 1121 |2E] for a single-phase membrane is defined as 



where H is the mean curvature of the membrane surface T, cq the spontaneous cur- 
vature, G the Gaussian curvature, a\ the surface tension, ai the bending rigidity and 
03 the stretching rigidity. 

In recent experimental studies, it has been found that the bending rigidity in the 
liquid-disordered phase differs from that in the liquid-ordered phase in two-component 
membranes [3 El- This can be attributed to, among other things, that the two phases 
have different lipid compositions or different concentrations of cholesterol molecules 
which serve as spacers between lipids. Thus, in the generalized bending elasticity 
model for two-component membranes, the bending rigidity is assumed to take 
values k\ and fc 2 respectively in two different components (phases) Ti C T and r 2 C T 
(with r = Ti U T 2 ). The common phase boundary between the two phases is denoted 
by 7o = Ti H r 2 . In general, the other parameters may also vary in different phases, 
however, in this paper, we ignore the effect due to a±, and Cq and concentrate only 
on the effect of the bending rigidity, though the methodology can be easily extended to 
more general case. In fact, the formulation we present here works for two components 
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having possibly different spontaneous curvatures, but for simplicity such curvatures 
are set to be zero in the numerical simulations. 

With the two phases co-existing on the membrane, it is natural to introduce a 
line tension on 70 to take into account the interfacial energy between the individual 
components [21 1221 I3T] . Coupling with the bending elastic energy, this leads to the 
following total energy determining the two component membrane 



where S is the line tension constant [221 • Note that in general, the line energy can 
also include the integral of a multiple of the curvature square on 70 [3] ■ 

The mathematical model that our study is based on is the minimization of the 
total energy defined in (|1.2() for a two component membrane with a prescribed to- 
tal volume, and prescribed surface areas of both components. In order to effectively 
model and simulate the experimental findings on the exotic morphology of the multi- 
component vesicles (mostly taken from [2]), we extend the recently developed phase 
field approach for the single component vesicles |18| to the multi-component case, 
which avoids the tracking of the vesicle membrane by viewing the surface and phase 
boundary as the zero level sets of phase field functions. The general phase field frame- 
work has been used successfully in many applications El El . For membrane de- 
formation, this approach has become increasing popular in the research community in 
recent years. So far, its applications have mostly confined to the case of using a single 
phase field function [51 ll5llTBll21H27j . albeit it is known that co-dimension two objects 
can be described effectively by a pair of level-set or phase field functions |?1 HtH l2l?l 132*] . 
With the introduction of a second phase field function, we demonstrate that the new 
two-component phase field model is capable of capturing rich complex morphological 
changes experimentally observed in the two-component vesicle membranes. Moreover, 
this model can be very easily generalized to study the open membranes or membranes 
with free boundary (see "3("j for experimental study and [3J "33J 1351 137] for analysis 
and computation). This is based on the observation that an open membrane can be 
thought as a two component vesicles with one component having zero bending rigidity. 
Further generalization is possible for vesicles with three or more components. 

The paper is organized as follows: in section 2, we present the phase field formu- 
lation of the total energy 1|1.2[1 and address the approach of penalty formulation for 
the constraints. In section 3, we first briefly discuss the discretization schemes and 
some implementation issue. After presenting some convergence tests to validate our 
method, we assemble a number of interesting experiments to explore the shape trans- 
formations due to the changes of different parameters. The numerical simulations are 
compared with experimental findings including the merge and splitting of different 
components. In section 5, we present the phase field formulation for open membranes 
and some numerical simulation results. We then make some conclusion remarks in 
section 6. Some technical derivations are provided in the appendix. 

2. A diffusive interface model. We start by introducing a pair of phase field 
functions (0(x), r;(x)), defined on the physical (computational) domain Q. 

The function (j> = <X X ) is used so that the level set {x : 0(x) = 0} gives the 
membrane T, while {x : </>(x) > 0} represents the interior of the membrane (denoted 
by fii) and {x : 4>(x) < 0} the exterior (denoted by Sl e ). In the phase field models of 
a single component vesicle, this is the only phase field function used |17|. 
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Next, we take another closed surface T± defined on domain Q and being perpen- 
dicular to r, such that it is the zero level set {x : 77 (x) = 0} of a phase field function 
r\ = r](x) in f2 with {x : r;(x) > 0} being the interior of T± and {x : </>(x) < 0} the 
exterior. We thus take the part of T in the interior of T± as the first component I\ 
and the remain part of T (denoted by 1^2) makes up the second component. Note 
that there may be many choices to select T±, but we are mostly interested in the 
level set {x : 7/(x) = </>(x) = 0} which gives the boundary between two components, 
with {x : T](x) > and 0(x) = 0} representing one component of the membrane and 
{x : r\ (x) < and <fi(x) — 0} the other component. 

In the phase field modelling, the functions r\ and <fi are forced to be nearly constant 
valued except in thin regions near the surfaces T and T± respectively. We use two 
small positive constant parameters e and £ to characterize the widths of the thin 
regions (also called the diffusive interfaces). We note that a phase field function 
(order parameter), like our r], has been introduced in |^ \T7\ to describe the phase 
segregation on the membranes, but different from our phase field description of the 
surface T, an explicit construction of the membrane surface and a direct computation 
of the bending elastic energy are used there instead of the phase field representation 
of the membrane surface. 

Similar to ^Hj, we have the phase field elastic bending energy defined by 

E(4>,v)= I ^{eA^+i-^ + coi^d-^Ydx, (2.1) 
isi 2e e 

where we take a variable bending rigidity given by k(rf) = k + ctanh(|), so that k + c 
corresponds to the value of the bending rigidity of one component and k — c the other. 
Similarly, 2co(?7) = ( c i + c 2) + ( c i — C2)tanh(|), so that c\ and C2 correspond to the 
spontanenous curvatures in the two components respectively. A few other functionals 
needed in our model are as follows: 

m V) = f n 4 |Vr?| 2 + ^(V 2 - 1) 2 ][||V0| 2 + \-y - \f\ dx, (2.2) 



A( ( f>)= J^m 2 + ^ 2 -l) 2 }dx, (2.3) 



A{(j>) = / 4>dx, (2.4) 
Jo. 

D&rj) = J tanh(|)[||V0| 2 + 1(0 2 - l) 2 ] dx . (2.5) 

To reveal the meaning of the above functionals, we follow similar discussions 
in ^S] to assume an ansatz of the form </>(x) ~ tanh(c?(x, T)/(y/2e)) and ?y(x) ~ 
tanh(c?(x, r^)/(v / 2C)) f° r the phase field functions. Here d denotes the signed distance 
function. In this ansatz, we can check that as e and £ tend to 0, that is, in the sharp 
interface limit, 

E{4>, rj) - ^ J2 1 k ^ H - c *) 2 ds ■ ( 2 - 6 ) 



More details are given in the appendix, along with a brief derivation of the folllowing 
asymptotic limits 

v{4>)^2\su\-\n\, A(<f>)^^\r\, (2.7) 

and 

L(4>,V) - \J Sdl , D(4>,ri) - - |r 2 |) . (2.8) 

To re-cap the discussion, the total energy in the phase field two-component model 

is 

£((f>,ri) = E(4>,ri) + Lfan) , 

while the constraints are given by 

V{4>) = v d , A{4>) = oo, D(<f>, 77) = a d , (2.10) 

with v d , at and a d being the prescribed volume difference (hence the interior volume is 
prescribed), the total surface area and the area difference between the two components 
(hence areas of both components are prescribed). 

To maintain the consistency of the phase field model which is based on <f> and 
77 having the tanh profiles and the orthogonality between T and T±, additional con- 
straints are imposed. First of all, the orthogonality constraint on the normal direc- 
tions of the two surfaces, written in our phase field formulations, can be enforced by 
V<p ■ V77 = on or near the phase boundary {x : <p(x) = jj(x) = 0}. With (f> and 77 
having tanh profiles, their gradients become small away from their zero level sets, the 
orthogonality constraint may thus be enforced everywhere by penalizing 

N(ct>,r]) = J^\V4>-V V \ 2 dx. (2.11) 

Secondly, to better maintain the tanh profile of 77, especially for the case with a large 
line tension energy, we have two options, one is to add a small regularization term, 
much like the bending elastic energy for <fi but with a very small bending rigidity; 
another option is to regularize through the following functional 

P(V) = f (f |Vr,| 2 - l(r, 2 - Iffdx , (2.12) 

which also vanishes for any function r\ with a tanh profile. 

Summarizing the above, the variational phase field model to describe the two- 
component vesicles in the energy minimizing state is to minimize the total energy 
£{<p, 77) = E(<f>, rj)+L((f>, 77) with constraints V{<p) = a\, A((f>) = a 2l D{4>, rj) = a 3 while 
N((f>,rf) and P(r?) remain small. So, by adding both the penalty and regularization 
terms, the vesicle surface and the two components phase boundary are determined by 
a pair of phase functions (</>, rj) which minimizes the energy 

£ M (4>, 77) = E(<f>, 77) + L(4>, 77) + \Mi{V{<l>) - v d f + \m 2 {A{4>) - a ) 2 

+±M 3 (D(<j>, V ) - a d f + ^M 4 (iV(^ r,)) 2 + ^M 5 (P(77)) 2 (2.13) 

where {Afj}f =1 are penalty constants for the constraints on the volume and surface 
areas while {Mi}f =4 are regularization constants for maintaining better control on the 
phase field functions. 
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3. Numerical Simulations of Two-Component Membranes. In this sec- 
tion, we compute the minimum of the phase field energy l|2.13|l by adopting a gradient 
flow approach which has been very effective for solving the phase field model of single 
component vesicles |17l 119) . The equations for the gradient flow are given by: 



The monotone decreasing of the energy Em is ensured for t > 0. For simplicity, we 
only consider the case where c\ = C2 = 0, this allows us to focus on examining how 
the variation in the bending rigidity alone affects the vesicle shape deformation and 
the equilibrium configurations of two-component membranes. The more general cases 
involving the spontaneous curvatures are to be considered in the future. 

Discretization and code development. . We take the spatial computational domain 
as the box £1 = [— tt,tt] 3 and assume that membranes are enclosed in the box. More- 
over, we choose to set £ — e in our numerical simulations. For the spatial discretization 
of H3.1(l in il, a Fourier spectral method is used. Due to the regularization effect of 
the finite transition layer, for fixed e and enough Fourier modes, the spectral method 
is an efficient way to solve (|3.1H with the help of FFT routines ^UJ- A couple of 
options are implemented for the time discretization, such as an explicit forward Euler 
scheme or a semi-implicit Euler scheme |17j . The time step At is chosen to ensure the 
energy decay. For most of our numerical experiments, although fully adjustable, At 
is kept in the range of 1CP 6 to 10~ 7 . The simulation codes are fully parallelized on 
both distributed memory systems via MPI and shared memory systems via OPENMP 
to improve its efficiency and functionality in conducting extensive three dimensional 
simulations. 

Problem set up and initial profiles.. We now discuss how we choose various pa- 
rameters in the simulations. Though in theory the gradient flow can be started from 
any pair of initial phase field functions, a proper choice often speeds up the evolution 
process and allows more efficient solution of the equilibrium state. With the penalty 
formulation, a particular constraint can be simply removed by setting the correspond- 
ing penalty constant zero. For example, setting M% = M3 = would eliminate the 
total area and area difference constraints. This fact can be utilized to find good initial 
phase field functions. 

For example, as illustrated in Figure lTTl for a given r > 0, we may start from 
two special phase field functions as </>(x) = tanh(-^jip) and rj(x) = tanh(-^) where 
z is the third component of x. This provides two hemispheres that represent the 
two components (colored in red and blue respectively, or in gray-scale represented 
by lighter and darker regions). Starting from this initial state, and setting Mi = 
to eliminate the volume constraint, the sphere gradually becomes more elliptical due 
to the presence of line tension, then further transform to a gourd like shape. We 
may stop at an intermediate shape and add back the volume constraint. This would 
provide a variety of initial shapes to be used in the simulations. 

Convergence verification.. For a particular numerical simulation, the quality of 
the numerical result may be affected by the choice of computational domain, the 
parameter e (the effective width of the diffusive interface), the number of grid points, 
and the choices of other parameters used in the simulation. The parameter e is 
generally taken to be a few percentage points of the domain size to ensure a relatively 
sharp interfacial region and the consistency with the sharp interface description (the 
e — > limit). The mesh size is normally taken to be several times smaller than the 



>t = - 




m = - 



Si] 



(3.1) 



5 



Fig. 3.1. Line tension drives a two- component sphere to a gourd shape. 
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Fig. 3.2. The 3d and cross-section views of tj> and the 3d view for part of n, computed with 
different parameters: left picture, e = 2h (above) and e = 1.5h (below); right picture, 64 3 grid 
(above) and 48 3 grid (below). 

width of the transition layer to ensure adequate spatial resolution. To ascertain the 
accuracy and robustness of our numerical algorithms and the parameter selections, 
we here present results of some numerical tests on the convergence and performance 
of our method. 

The first set of experiments given in Figure is designed to test the dependence 
of the resolution of the phase field function on the the parameter e and grid size. We 
take a shape similar to the previous experiment. First, we take a 64 3 grid but use 
different values of e at 0.1964(= 2h) and 0.1472(= 1.5h). The other parameters are 
defined by v d = -216.52, a = 29.46, a d = 0.23 and M* = 3.2 x 10 5 for all i. The 
two equilibrium shapes are almost the same except the transition layer width. The 
corresponding final energy values 124.49 and 123.82 are very closed to each other. 
The left picture of Figure 13.21 gives the final three dimensional views and some cross 
section views of the phase field functions <f> and r\. 

Now we use the same set of parameters (e = 0.1964, same initial <po in the same 
domain), but solve the problem on two different grid sizes 48 3 and 64 3 . We set the 
parameters v d — —216.52, ao = 29.46, a d = 0.23 and constants Mj = 10 4 for all i. 
The right picture of Figure 13.21 provides the details of the simulations, with the 3d 
views of <f> and their density plots of the cross-sections in x — z plan. The plots of the 
corresponding rj are similar to that in the third column of the left picture of Figure 
13.21 and are thus omitted. The final values of energy are 124.39 and 124.42 while 
the elastic bending energy values are at 48.05 and 47.96, and the line tension energy 
values at 76.34 and 76.46 respectively. The close values substantiate the convergence 
of the simulated results. 

The convergence can also be verified for different penalty and regularization con- 
stants {Mj}| =1 . The difference in the penalty and regularization is to be understood 
as follows: the penalty constants {Mi}^ =1 are taken to be larger and larger to ascer- 
tain the satisfaction of the volume and areas constraints. The regularization constants 
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Convergence of the Lagrange multipliers. 
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Table 3.2 

The diminishing effect of regularization on the total energy. 



{M4, M5}, on the other hand, are taken to be smaller and smaller so that while the 
orthogonality of the zero level sets of the two phase field functions and the tanh like 
profile of r\ are both effectively maintained in the simulations, the associated energy 
contributions from the regularization terms in fact diminish. 

First, we define the Lagrange multipliers by Xi = limM^oo A,({Mj}f) with Ai = 
Mi(V(4>)-Vd), A 2 = M 2 (A(4>)-a ), A 3 = M 3 {D{<f>,>r]) - a d ). With other parameters 
given {Mi = M 5 = 10000, v d = -216.52, a = 29.46, a d = 0.230, e = 1.768, 
h = 0.17355), we set larger and larger values for M\ = M 2 = A/3. The results 
are given in Table ETH which show that Ai, A 2 and A 3 converge to the Lagrange 
multipliers, and errors in constraints also decrease. 

Next, we demonstrate that the regularization terms provide effective control on 
the phase field functions but do not contribute significantly to the energy minimiza- 
tion. We set a sequence of decreasing values for M4, M5 while taking the same values 
for Mi = M 2 = M3 = 10000, and keeping the values of other parameters the same as 
in the previous test. The results are given in Table l3~2l where £4 = ^Mi{N{4>, rj)) 2 
and E 5 = ^M 5 {P{r/)) 2 and their ratios with the total energy £m are provided. We 
can observe the diminishing and negligible effect of the regularization terms while 
there is no noticeable change in the simulated membrane. 

Having demonstrated the convergence of the numerical algorithms, we next study 
the effect of different bending rigidities and various line tension constants. Then by 
adjusting the bending rigidities in the two components and the line tension, we can 
simulate the the vesicle shapes in experiment findings [2]. Unless noted otherwise, 
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the simulation results reported in the following are obtained with a 64 3 grid sizes and 
e = 0.1736 which can provide sufficient resolution based on the convergence study. 

Effect of the bending rigidities. . The values of bending rigidities often play a key 
role in forming various shapes of vesicles. Our first experiment is a simulation of the 
striped vesicles. We start from an initial shape where the red component is situated 
in the center to give a stripe-looking vesicle. As shown in the first row of the Figure 
13.31 with parameters Vd = —213.98, ao = 29.46 and a^ = —13.31, the initial shape 
grows into a very regular stripe-looking ellipsoid shown in the middle of the first row. 
In this experiment, the bending rigidity for the red component is 1.0 whereas the blue 
component is 3.0. With line tension being fixed at 10.0, we then make a switch of the 
bending rigidity of the two components. As shown in the right picture of the first row, 
the red component of the ellipsoid in the middle grows to a more cylindrical like shape. 
Next, by preserving the bending rigidity of the blue component while increasing that 
of the red component from 3.0 to 19.0, the red component shrinks in the middle 
and we get a thinner center band as shown in the left picture of the second row. It 
is obvious that the concave region has a smaller mean curvature. We can further 
increase the area of the blue component by setting ao = 33.46 and a^ = —19.82, and 
with bending rigidities 3.8 and 0.2 respectively for the red and blue components, we 
get the middle picture of the second row in Figure 13.31 One can compare it with 
the last picture found in actual experiments [2] though the differences of the bending 
rigidities are not as significant as those used here. In the final shape, the center band 
has the lowest mean curvature and it is occupied by the red component (having larger 
bending rigidities). 




Fig. 3.3. Different values of bending rigidity lead to different shapes of striped vesicles (the 
right bottom picture is reproduced from 

As expected, the numerical simulation shows that the component with a larger 
bending rigidity is more likely to remain in regions with smaller values of mean cur- 
vature. 

Effect of line tension constants.. By intuition, we expect that larger line tension 
generally leads to a shorter interfacial line between two different components. And the 
line tension is balanced by the bending and elasticity force and the volume constraint. 
In most of the cases, the volume constraint plays a key role in balancing a large line 
tension as in the experiments illustrated by Figure l3~4l and l3~Bl 

In Figure 13.41 the pictures shown there correspond to equilibrium shapes with 
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Fig. 3.4. Different values of line tension result different vesicle shapes. 

three different values of the line tension 10.0,30.0,100.0. The bending rigidities of 
the blue colored component is 3.5 while that of the red is 0.5. By increasing the 
line tension, the individual components in the two-component vesicle become more 
hemisphere like which are the results of the increasing effect of line tension under the 
same volume and surface area constraints. 

Figure 13.51 gives an even more convincing example to the rupture and vesicle 
fission observed in this process. As shown in Figure 1331 we start from the top left 
shape. While preserving Vd, a-o and ad to be —213.98, 29.46 and —13.31 respectively, 
we increase significantly the line tension from 10.0 to 100.0. The vesicle gradually 
breaks its vertical symmetry and a small blue vesicle is separated and eventually 
absorbed into the top portion through a process like Oswald ripening. Finally, the 
vesicle (bottom-right picture of Figure I3.5fl only contains two parts, much like the 
shape observed in the experiments 0j. 



• P 



t ' 




Fig. 3.5. Effect of line tension: rupture and fission of vesicles components (the right bottom 
picture is reproduced from Wl). 



Comparison with other experimental results.. We now focus on the simulations 
that mimic other two-component vesicle shapes observed in the experiments of 
similar to the results depicted in Figures 1331 and 1331 

As shown in the two rows of Figure l33l we carry out two simulations starting from 
a shape given on the left. In both simulations, the red component has bending rigidity 
3.0, and the blue component has bending rigidity 1.0. The line tension between two 
components is 30.0. The parameter Vd for volume constant is —218.00, and the surface 
area parameter ao is 29.46. The parameter ad giving the difference of surface areas 
of the two components takes on the values 18.76 and —18.76 respectively. The final 
shapes of the two simulations are shown in the center pictures in both rows. One can 
compare them with the right most experimental picture provided in 
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Fig. 3.6. Similar membrane shapes with different areas for the two components (the pictures 
on the right column are reproduced from H?T). 
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Table 3.3 

Energy comparison for the shapes given in Figure T3. 61 



The energy values of the two experiments illustrated in Figurc EQSl are given in Ta- 
blc !3.3l with energy contributions listed for individual components and the line tension 
from the phase boundary. We get almost the same line tension energy contribution, 
but, as caused by the difference in the bending rigidities, the elastic bending energy 
contributions differ by a factor of 3, which is reflective of the ratio of the bending 
rigidities. 

We now turn to simulate a couple of other interesting shapes experimentally 
observed in as illustrated in the last pictures of Figure 13.71 and Figure 13.101 re- 
spectively In Figure 13.71 we first start from a spherical surface which is divided 
into two components where one component occupies similar spherical caps in twelve 
well-spaced locations on the membrane surface. The bending rigidity is 3.5 for the 
red component and 0.5 for the blue component, and the line tension is 10.0. With a 
larger surface area of the blue component and a smaller volume than those values for 
the exact sphere, the blue component (with smaller bending rigidity) starts to bulge. 
The parameters are taken respectively as vj, — —174.17, clq = 54.63 and a ( i — 11.01. 
If we further increase the volume and enlarge the relative area of the blue compo- 
nent by increasing Vd to —167.00, while keeping ao at 54.63 and changing a<j to 5.00, 
the resulting computed shape (the third picture in Figure l3~T|) is very similar to the 
experiment findings |2j (the last picture in Figure l3~7jl . 

The shape corresponding to the third picture of Figure l3~7l stavs as a near equilib- 
rium (meta-stable) state for a range of parameter values. But if we further increase the 
area of the blue component, for example, by setting a c i = 2.5, further coarsening of the 
blue components will take place. The merger of disconnected components continues, 
much like the Oswald ripening effect, and eventually transforms into shapes similar 
to that presented earlier in Figures 13.21 and 13.41 The transformation is illustrated in 
Figure El 
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Fig. 3.7. A sphere with disk like bumps: comparing with biological experiments. 
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Fig. 3.8. The merging of blue component (cut view). 



Next, we take an initial membrane profile similar to that in the second picture of 
Figure ET7I By setting Vd = —203.00, a = 54.63 and ad — 5.00 so that both the total 
volume and the area of the red components are decreased, we can then observe the 
growth of bumps of the red component, leading to a shape shown in the right pictures 
of Figure Take other initial profiles, other equilibrium shapes as shown in the left 
and center pictures in Figure 1531 have also been observed in our simulations. 




Fig. 3.9. Various shapes of two component membranes. 



From FigureEnH it can be seen that the two-component vesicles may display very 
rich patterns, even in the absence of spontaneous curvature effect. One naturally may 
wonder if some of them are experimentally observable. The next set of experiments 
draws inspiration from the center and right figures of Figure 13.91 and leads to inter- 
esting comparisons with similar experimental observations in ;2J. We start with the 
same phase field <fi as the profile in the right picture of Figure l3"31 but use a modified ij 
such that the neck of the bumps are formed by the blue component as the case of the 
center picture of Figure l531 This leads to an initial shape as shown in the left picture 
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of Figure IS. 1UI Setting the parameters as — —203.00, a = 54.63, and = 22.68, 
we finally get a shape (center picture of Figure lS.lOfl very close to the experimentally 
observed shape given in |2] (right picture of Figure I3.10fl . 




Fig. 3.10. Two-component shape with 14 bumps (the last picture is from T§f ). 

Shapes depicted in !3.10l are fairly robust. In fact, with a slight modification of the 
final shape and a rotation with a given angle, then we found that using the gradient 
flow, the equilibrium solution is again in the shape (except for a rotation). Results of 
such calculations on both 64 3 and 96 3 grids are given for comparison. 




Fig. 3.11. Rotated two-component shapes with 14 bumps computed by different meshes. 

4. Open liposomal membranes. In this section, we apply similar ideas to 
model open lipid membranes. The transformations from vesicles to open membranes 
and the reverse process from open membranes to vesicles were first observed in |30| . 
Here, we only consider the one-component open membranes with specified surface 
areas. The total energy of an open membrane T with edge 70 may be conveniently 
defined as the sum of the elastic bending energy and the line tension energy [HI 1331 

EHESIEZJ: 

/ (ai + a 2 (H - c ) 2 + a 3 G) ds + / Sdl. 
Jr J l0 

For simplicity, we set the surface tension ai and the line tension 5, as two constants, 
we also do not consider the contribution of the geodesic curvature term in the line 
energy on the boundary. The effects of surface tension, the Gaussian and spontaneous 
curvatures are also ignored. Our problem is then to minimize the following total 
energy 

E = kH 2 ds+ / Sdl 
Jr J jo 

with prescribed surface area |T|. 

Most of the available numerical simulations for open membranes have largely 
been confined to axis-symmetric cases based on the variational calculation of the 
above energy. We hereby develop a new phase field model for open membranes, and 
present some numerical simulations for the full three dimensional case to demonstrate 
the effectiveness of the model. 
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Phase field model for open membranes.. We can treat open membranes as two 
component membranes with one component having zero bending rigidity. Again, we 
let 70 be the intersection of two orthogonal surfaces Y and Tq which are implicitly 
defined as the level-set of the functions <j> and r\ respectively. 

Now we denote c(rf) = |(1 + tanh(|)), and let the line tension energy be still 
formulated by L((f>, rj) in (|2.2p . with the elastic bending energy of the membrane 

Then, our phase field model for open membranes is to minimize E((f>, rj) +L((j), rf) with 
the surface area constraint 

D(4>,r,)= J c(r,)[^\Vcf>\ 2 + ^ 2 -l) 2 }dx = a d . (4.1) 

Similar to the two-component vesicle case studied earlier, to maintain the good profiles 
for both phase field functions 4> and n and the orthogonality of T and T± , we can again 
take the penalty formulation 

£ M {4>, V) = WW, r,) + L{ct>, V ) + \m^D{4>, V ) - a d f 

+ l -NU{N{^ 1 )) 2 + ~M 5 (P( V )) 2 + lM 6 (P(<b)) 2 . (4.2) 

Then, we can again use a gradient flow like l|3.1|l to compute the equilibrium shapes 
by a similar numerical scheme as that given in section [3] 

Numerical simulations of open membranes. We now present some numerical sim- 
ulations of open membranes and compare them with biological experimental findings. 
Most of the model and simulation parameters are chosen to be in the same range as 
that for the two component vesicle simulations in the earlier section. 

83 

Fig. 4.1. Open membranes with different line tensions (the right most picture is reproduced 
from HUf ). 

Figure 14.11 gives the simulation results of a simple open membrane. Starting 
from a half sphere (the left picture), with bending rigidity k = 1.0 and line tension 
S = 1.0, we get an equilibrium shape shown in the second picture. If a larger line 
tension 8 = 1.28 is used, an equilibrium shape is reached as that in the third picture. 
One can compare it with the right picture obtained in the biological experiments 
described in (3U|. We note that the elastic bending energy are 10.12 and 15.96 and 
the line tension energy are 10.94 and 7.95 respectively for the solutions in the second 
and third pictures. 

The time evolution snapshots are given in Figure 14.21 where the line tension is 
taken as S — 25.0. The simulation results show that, when the line tension becomes 
large enough, the open membrane becomes self-enclosed. 
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Fig. 4.2. Open membrane closes due to large line tension. 



Finally, in Figure EHfl we simulate a shape (the right picture) with three holes as 
observed in an experiment of |3U| . Starting from the left most picture corresponding 
to an ellipsoid with three holes, setting the bending rigidity k = 1.0 and line tension 
5 = 1.0, and following the gradient flow of the energy, the initial shape starts to 
deform first into an intermediate shape given in the second picture. The computed 
equilibrium shape is shown in the third picture which again shows striking similar to 
the experimental finding. 




Fig. 4.3. Open membranes with three holes (the right most picture is reproduced from \3Ut ). 

5. Conclusion. In this paper, we formulated a phase field model for the multi- 
component vesicles membranes, and as a special case, the open membranes with free 
edges. The models incorporate the effect of the elastic bending energy together with 
the line tension between each two components. Full three dimensional numerical 
simulations presented here demonstrate that the experimental observations given in 
[2] can be effectively simulated by the phase field bending elasticity and line ten- 
sion model. Furthermore, the simulation results illustrate that many experimentally 
observed exotic patterns such as bud formation and vesicle fission can appear in two- 
component vesicles due to the inhomogeneous bending stiffness and the competition 
of the bending energy and the interfacial line tension even without incorporating the 
spontaneous curvature or the asymmetry of the bilayer. 

In conclusion, we point out the this generalization of our diffusive interface model 
to two component vesicle membranes fits nicely into the previously established unified 
framework for the derivation of dynamic and static equations and the development 
of numerical algorithms and codes. Many issues remain to be examined in future 
works. First, there maybe other more effective ways in formulating the line tension 
energy, including the use of geodesic curvature along the phase boundary, which may 
be important to model the difference of the stretching rigidity in two components 
[2| • Second, more rigorous analysis of our models are needed in the future. Third, 
in our numerical simulations, we have not examined the effect of the spontaneous 
curvature as we have done for the one component case |15| . It is expected that more 
complex shapes would be be discovered in this case. Moreover, the interaction of 
multi-component vesicles with the fluid and electric fields are also exciting topics to 
be investigated further in the future. 
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Appendix: Justification of the energy and constraints.. We now provide 
some brief calculations to rationalize the definitions of the energy functional and the 
constraints in the phase field setting. Same as the discussion in JT7], we can first 
illustrate that in a very broad ansatz, for small e and £, minimizing W(<j), rf) leads to 
a phase field function <j>(x) which is approaching to tanh(d(a;, T)/(^/2ej) as e — > 0. In 
fact for small e, due to the uniform bound of the functional B, the region away from 
the level set <fi = are all close to <f> = +1 or —1. In such cases, one can always define 
the following transformation near the interface: 



0) = q (^)i 



(5.1) 



where d(x) is the distance of the point x £ O to the surface T. Substituting this into 
Q2.1[l . we have that: 



E{4>) = 



2e 



q"(^-)Ad(x) + k<l £ "-(<l e2 -l)<l e ) 
e e 



dx. 



(5.2) 



If we keep k(rf) positive, as e — > 0, to minimize the energy, the leading term in the 
above has to vanish, that is, 



ell (e2 



(<T - i)g e 







(5.3) 



which means that the transition region profile q e (-) is approaching to the function 
tanh(-^). In the meantime, we see that <j) is approaching to the Heaviside function 
with 1 inside of the interface and —1 outside. V still coincides with the zero level set 
of (p. Moreover (|5.1I) indicates that the parameter e is effectively the thickness of the 
transition region between {</> = 1} and {<fi = —1}. One can refer |14j for more rigorous 
proof of this. 

Now, we denote s((j>) = ||V</>| 2 + j^{<t> 2 — l) 2 - When the line tension L(<fi, rj) reach 
its minimum, we have 



SL 
Si] 



1 



fir], s) = -eV ■ (sVr?) + -s{if - 1)?/ = 



For <j) = tanh(d/(\/2e)), s = ^1 
Then /(rj, s) = means 



l) 2 and therefore Vs • Vrj = as V0 • V?/ = 0. 



-eArj + - (rf - l)r) = . 
If we write rj again by q e (d(x ,T j_) / 'e) , from the above equation we have 



- q ^(-)Ad(x) + -((q^-l)^ 



As e — > 0, we have (q t2 — l)q e — q fJI — 0. To minimize L, we can expect that far 
away from the T±, q e is +1 or — 1, therefore we also have q £ (x) — tanh(^^). On 
the other hand, we can use the same argument for (f> if we know r\ is a tanh function, 
which would further strengthen the ansatz that <p and rj are both tanh functions to 
lead order of e. In fact, following more careful analysis as those in we expect 
that the differences between <p and rj and the respective tanh profiles are second order 
in e which would allow us to rigorous derive the asymptotic limits (|2.6I2.8|) . 
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